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Abstract 

The partition function of a quantum field theory with an exact symmetry can be de- 
composed into a sum of functional integrals each giving the contribution from states 
with definite symmetry properties. The composition rules of the corresponding transfer 
matrix elements can be exploited to devise a multi-level Monte Carlo integration scheme 
for computing correlation functions whose numerical cost, at a fixed precision and at 
asymptotically large times, increases power-like with the time extent of the lattice. As a 
result the numerical effort is exponentially reduced with respect to the standard Monte 
Carlo procedure. We test this strategy in the SU(3) Yang-Mills theory by evaluating 
the relative contribution to the partition function of the parity odd states. 



1 Introduction 



Dynamical properties of quantum field theories can be determined on the lattice by 

computing appropriate functional integrals via Monte Carlo simulations. For the most 
interesting theories this is, up to now, the only tool to carry out non-perturbative 
computations from first principles. The mass of the lightest asymptotic state with 
a given set of quantum numbers can, for instance, be extracted from the Euclidean 
time dependence of a suitable two-point correlation function. Its contribution can be 
disentangled from those of other states by inserting the source fields at large-enough 
time distances. The associated statistical error can be estimated from the spectral 
properties of the theory [1,2]. Very often the latter grows exponentially with the time 
separation, and in practice it is not possible to find a window where statistical and 
systematic errors are both under control. This is a well known major limiting factor 
in many numerical computations such as, for example, the computation of the glucball 
masses in the Yang-Mills theory. A widely used strategy to mitigate this problem is to 
reduce the systematic error by constructing interpolating operators with a small overlap 
on the excited states [3,4]. The lowest energy is then extracted at short time-distances 
by assuming a negligible contamination from excited states, sometimes also with the 
help of anisotropic lattices [5,6]. This procedure is not entirely satisfactory from a 
conceptual and a practical point of view. The exponential problem remains unsolved, 
and the functional form of the sources are usually optimized so that the correlator shows 
a single exponential decay in the short time range allowed by the statistical noise. A solid 
evidence that a single state dominates the correlation function, i.e. a long exponential 
decay over many orders of magnitude, is thus missing. 

In this paper we propose a computational strategy to solve the exponential prob- 
lem. The latter arises in the standard procedure since for any given gauge configuration 
all asymptotic states of the theory are allowed to propagate in the time direction, re- 
gardless of the quantum numbers of the source fields. By using the transfer matrix 
formalism, we introduce projectors in the path integral which, configuration by config- 
uration, permit the propagation in time of states with a given set of quantum numbers 
only. The composition properties of the projectors can then be exploited to implement 
a hierarchical multi-level integration procedure similar to those proposed in Refs. [7,8] 
for the Polyakov loops. By iterating over several levels the numerical cost of computing 
the relevant observables grows, at asymptotically large times, with a power of the time 
extent of the lattice. 

We test our strategy of a "symmetry constrained" Monte Carlo in the SU(3) Yang- 
Mills theory by determining the relative contribution to the partition function of the 
parity-odd states on lattices with a spacing of roughly 0.17 fm, spatial volumes up 
to 2.5 fm^, and time extent up to 3.4 fm. The algorithm behaves as expected, and in 
particular the multi-level integration scheme achieves an exponential reduction of the 
numerical effort. In the specific numerical implementation adopted here the computa- 
tion of the projectors is the most expensive part, and its cost scales roughly with the 
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square of the three-dimensional volume. The realistic lattices considered in this paper, 
however, were simulated with a modest computational effort. 

The strategy proposed here is rather general and we expect it to be applicable to 
other symmetries and other field theories including those having fermions as fundamen- 
tal degrees of freedom. It can, of course, be quite useful also for computing excited 
levels in other quantum mechanical systems. The basic ideas were indeed checked in a 
considerable simpler and solvable quantum system with a non-trivial parity symmetry, 
namely the one dimensional harmonic oscillator [9]. 



2 Preliminaries and basic notation 



We set up the SU(3) Yang-Mills theory on a finite four-dimensional lattice of volume 
V = T X with a spacing a and periodic boundary condition^. The gluons are 
discretized through the standard Wilson plaquette action 



S[U] 



^^[l-^ReTr{[/^,(x)} 



(2.1) 



where the trace is over the color index, /3 = G/g^ with qq the bare coupling constant, 
and the plaquette is defined as a function of the gauge links U^{x) as 



Uf,uix) = Uf,{x) Uy{x + A) UUx + u) Ul{x) , 



(2.2) 



with /i, z^ = 0,...,3, /iis the unit vector along the direction ^ and x is the space-time 
coordinate. The action is invariant under a gauge transformation 



U^{x) U^{x) = n{x) U^{x) ^\x + (i) 
with Q.{x) G SU(3). The path integral is defined as usual 



(2.3) 



Z = / Di[U] e 



-S[U] 



X i_i=0 



(2.4) 



where T>U is the invariant Haar measure on the SU(3) group, which throughout the 
paper will be always normalized such that / DU = 1. The average value of a generic 
operator O can thus be written as 



(O) 



1 

z 



(2.5) 



^Throughout the paper dimensionful quantities are always expressed in units of a. 



2 



2.1 Hilbert space 



The Hilbert space of the theory is the space of all square- integrable complex- valued func- 
tions tplV] of Vfc(x) G SU(3) with a scalar product defined as (x is the three dimensional 
space-coordinate and k = 1,2,3) 

3 

(cl^li;) = / J^s[V]^[V]*^[V] , D3[y] = n U^Vki^) ■ (2-6) 

The "coordinate" basis is the set of vectors which diagonalize the field operator at all 
points X, i.e. 

Vfc(x)|y) = Vkix)\V) , (2.7) 



and which are normalized such that 



(2.8) 



From a quantum mechanical point of view, the field values (x) form the set of quantum 
numbers that label the vectors of the basis. In a gauge theory physical states are wave 
functions which satisfy 

i^iV"] = (2.9) 
for all gauge transformations 0. A projector onto this subspace can be defined as 

(F|Pg|V> = / D[n] V'lO , D[0] = llvn{^) , (2.10) 

X 

and it is straightforward to verify that Pq = Pq. 
2.2 Transfer matrix 

The transfer matrix of a Yang-Mills theory discretized by the Wilson action has been 
constructed many years ago [10-13]. The subject is well known and it appears on text 
books, therefore we report only those formute which are relevant to the paper. The 
starting point is to rewrite the functional integral in Eq. (|2.4p as 



T—l 

Z = [ U D3[FJT 
where the transfer matrix elements are defined as 



with 



















L 


VXQ + I, Vxg 


= K 


Kco+l) Vxo 


+ \w 






Vx,_ 



(2.11) 

(2.12) 
(2.13) 
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and O being identified with the hnk in the temporal direction. The kinetic and the 
potential contributions to the Lagrangian are given by 



K 



and 



w 



^Yl l-^ReTv{yfc(xo + l,x)yJ"(xo,x)} 

x,k 



XO 



l-^ReTi{l4,(io,x)} 



X k,l 

respectively, where V^i is the plaquette defined in Eq. 



(2.14) 



(2.15) 



computed with the links 



XO 



w 



"xo 



, while the depen- 



Vfc(x). The potential term is gauge- invariant, i.e. W 

dence of the kinetic term on the gauge transformations at time (xq + I) and at time 
xq is only via the product 17^0'. Thanks to the invariance of the Haar measure under 
left and right multiplication, this implies that the transfer matrix is gauge-invariant 



yQ' yQ 
•^Xo + l' XO 



X0 + I7 ^Xo 



and that 



(2.16) 



(2.17) 



The latter are thus matrix elements of a transfer operator T between gauge invariant 
states 



T 



Vxo + llVxg 



F^o+iIPgTPgIF, 



XO / ! 



and the functional integral can then be written as 

. 1 T 



Tr 



TP 



G 



(2.18) 



(2.19) 



where the trace is over all gauge invariant states. For a thick time-slice, i.e. the ensemble 
of points in the sub-lattice with time coordinates in a given interval [xq, uq] and bounded 
by the equal-time hyper-planes at times xq and yo, the transfer matrix elements can be 
introduced by the formula 



Vyg , VxQ 



2/0-1 



yo-i- 

n D3[^-o 

uio=xo+l zo=xo 



n -^[y^ 



(2.20) 



3 Decomposition of the functional integral 

The invariance of the system under a global symmetry can be exploited to decompose 
the partition function into a sum of functional integrals each giving the contribution 
from states with definite symmetry properties. In the following we will focus on the 
invariance of the Yang-Mills theory under parity. 
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In the coordinate basis, the parity transformation on gauge invariant states can be 
defined as 



P |v) = |v^) , |v) = Pg|f) , If (x) = - k) , 

which imphes that p =1. The parity eigenstates can then be written as 



|V) ± |V^ 



p|V,±) = ±|V,±) 



and their transfer matrix elements are given by 



T 



^0 + 1 ' ^0 



+ sT 



The invariance of the action yields 



T 



T 



and therefore 



Vxo+l,V,^o 



V"^ V 



Vxo+uV^^o 



Vxo+i,VX, 



(3.1) 

(3.2) 

(3.3) 
(3.4) 

(3.5) 
(3.6) 



For a thick time-slice the matrix elements between parity states can be introduced by 
exploiting the composition rule 



V V 

where xq < ^^o < Ho and in general 



T* 



5 ^0 



r-£S 



Vz,,Vx, 



} = j D3[y.o 



^0' ^0 



Vz,,Vx, 



(3.7) 



(3.8) 



It is easy to show that, in addition to relations analogous to those in Eqs. (j3.4p ~ (j3.6p . 
the identities 



T" 



^0 5 ^0 



hold. In particular they imply that 



T 







^0 5 ^0 



2:0 J 



D4[C/]sub e" 



-S[C/] 



T^[C/,„,t/,„_ij 



(3.9) 
(3.10) 

(3.11) 



an useful expression for the practical implementation of the multi-level algorithm de- 
scribed in the following section. The subscript "sub" indicates that the integral is 
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performed over the dynamical field variables in the thick time-slice [xo,yo] with the 
spatial components Uk{x) of the boundary fields fixed to Vfc(xo,x) and Vk{yo,5l) respec- 
tively. Finally, by inserting Eq. (|3.4p into Eq. (|2.1ip and repeatedly applying Eq. (13. 9p . 
it is possible to rewrite the path integral as a sum of functional integrals 



s=± 



T-1 

Z'= I H D3[Ko]T^[^xo+i,K.o 
xo=0 



(3.12) 



each giving the contribution from gauge-invariant parity-even and -odd states respec- 
tively 



-EoT 



n=l 



z- 



-EoT 



m=l 



-Err,T 



(3.13) 



In these expressions Eq is the vacuum energy, and E^ are the energies (with respect 
to the vacuum one) of the parity even and odd eigenstates, and and are the 
corresponding weights. The latter are integers and positive since for the Wilson action 
the transfer operator T is self-adjoint and strictly positive [11]. 

It is interesting to notice that even though the transfer matrix formalism inspired 
the construction, the above considerations hold independently of the existence of a 
positive self-adjoint transfer operator. The insertion of T*[V^(,, Vxq\ in the path integral 
plays the role of a projector, as on each configuration it allows the propagation in the 
time direction of states with parity s only. Indeed the parity transformation of one 
of the boundary fields in T[yj^(,, V^q] flips the sign of all contributions that it receives 
from the parity-odd states while leaving invariant the rest. The very same applies to 
the path integral in Eq. (j2.4p if the periodic boundary conditions are replaced by p- 
periodic boundary conditions, i.e. Vt = Vq . All contributions from the parity odd 
states are then multiplied by a minus sign. Similar considerations have already been 
exploited in different contexts, for instance in the study of the interface free energy of 
the three-dimensional Ising model [14]. 



4 Multi-level simulation algorithm 

The composition rules in Eqs. ()3.7p - ()3.10p are at the basis of our strategy for computing 
Z^ jZ (as well as a generic correlation function) with a hierarchical multi- level integration 
procedure. 



4.1 Projector computation 

To determine the parity projector between two boundary fields of a thick time-slice, the 
basic building block to be computed is the ratio of transfer matrix elements 
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The parity transformation in the numerator changes one of the boundary fields over the 
entire spatial volume of the corresponding time-slice, a global operation which could 
make the logarithm of this ratio proportional to the spatial volume, see for instance [14]. 
The transfer matrix formalism and the expected spectral properties of the Yang-Mills 
theory however suggest that, in a finite volume and for d large enough, only a few of 
the physical states give a sizeable contribution to this ratio, which is therefore expected 
to be of 0(1). These general properties can be studied analytically for the free lattice 
scalar theory, see for instance [15]. It goes without saying that the latter has a different 
spectrum from the Yang-Mills theory, and therefore can be used only as an example 
where our strategy can be studied analytically. 

Even tough the ratio R is expected to be of 0(1), the integrands in the numerator 
and in the denominator on the r.h.s of Eq. (|4.ip are, in general, very different and the 
main contributions to their integrals come from different regions of the phase space. 
The most straightforward way for computing R is to define a set of n systems with 
partition functions Zi . . . Z^ designed in such a way that the relevant phase spaces of 



successive integrals overlap and that Z\ = T[ya.(,+(i, Vjp] and Z^ 
ratio R can then be calculated as 



Z\ Z2 
R- = ^ X ^ X 
Z2 Zj, 



X X 



Z, 



n-1 



Zr, 



T[y,o+rf>^-o]- The 



(4.2) 



with each ratio on the r.h.s. being computable in a single Monte Carlo simulation 
by averaging the proper reweighting factor. To implement this procedure we start by 
generalizing the definition of the transfer matrix element in Eq. (j2.17p as 



D[1^']D[17] e 



XQ+l 



(4.3) 



where r G [-1/2,1/2] and 



+ r] K 



+ { --r]K 



V 









Vxo + l 







+ -w 



Analogously, Eq. (IS.llh can be generalized as 



(4.4) 



XQ+d—l 

n D3[K,o] 

W0=X0 + 1 



'xo+d-2 

n T 

. Z0=X0 



T 



yxo+d,Vxa+d-i,r 

(4.5) 



and the ratio R[T4Q+(i, Vx^] can be written as 

R[Vx,+d, Vx,] = n Wxo+d, ^xo, -1/2 + {k- 1/2) e] 



(4.6) 



k=l 
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where _ 

[ xo, J T[Ko+d,'^xo,r+e/2] ^ ' 

and e = With this choice of e the relevant phase spaces of two consecutive 

integrals overlap since the actions differ by a quantity of 0(1), while their fluctuations 
are of 0{^/V). To compute each ratio on the r.h.s. of Eq. (j4.6p one starts by noticing 
that the group integrals on Q' and Q. in Eq. (|4.3p can be factorized by introducing on each 
point of the time-slice xq the usual temporal link Uq{xq,x) = Q){x)VL (x) and a second 
temporal link U4^{xq,x) = Vl\x)^ {~^)- The average of the reweighting factor is then 
computed with the three-level algorithm described in Appendix[Al As other methods for 
computing ratios of partition functions which are present in the literature [16-18], the 
numerical cost scales roughly quadratically with the three-dimensional volume. Since 
the main goal of this paper is to present and test the validity of the strategy, we leave 
to future studies the development of a more refined and better scaling algorithm for the 
computation of the projector. 

4.2 Hierarchical integration 

Once the projectors have been computed, the ratio of partition functions /Z can be 
calculated by implementing the hierarchical two-level integration formula 



where ^ vq^xq 



Z' _ 1 
is defined as 



m,d 



D4[C/]e-^[^lp^^,ir,0l (4.8) 



T-T T'^[^3;o + (t+l)-d; UxQ+i-d] 
"^[Uxo + ii+iydiUxo+i-d] 



i=0 



with m > 1 and yo = XQ + m-d. The procedure can, of course, be generalized to a multi- 
level algorithm. For a three-level one, for instance, each ratio on the r.h.s of Eq. (j4.9p can 
be computed by a two-level scheme. Thanks to the composition rules in Eqs. (|3.7p and 
(j3.9p . the r.h.s. of Eq. (j4.8p does not depend on m and d. When computed by a Monte 
Carlo procedure, however, its statistical error depends strongly on the specific form of 
d 2/0) 3^0 chosen. The algorithm therefore requires an optimization which in general 
depends on the spectral properties of the theory. It is however important to stress that 
the multi-level hierarchical integration gives always the correct result independently on 
the details of its implementation. This can be shown by following the same steps in the 
Appendix A of Ref. [8]. There are two main differences: auxiliary link variables and 
their own actions need to be introduced for each value of r, and the computation of R 
requires a thermalization procedure for each value of r. We do not expect the latter to be 
particularly problematic since, as mentioned earlier, expectation values for consecutive 
values of r refer to path integrals with the relevant phase spaces which overlap. The 
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ratios R are computed by simulating systems corresponding to consecutive values of r 
one after the other, and by starting from the one used to extract the boundary fields 
(r = 0.5). 

4.3 Exponential error reduction 

The statistical variance of the estimate of a two-point correlation function (O(xo)O(O)) 
of a parity-odd interpolating operator O, computed by the standard Monte Carlo pro- 
cedure, is defined as 

= {0\x,)0\Q)) - (0(xo)0(0))2 . (4.10) 

At asymptotically-large time separations the signal-to-noise ratio can be easily computed 
via the transfer matrix formalism which, for <C xq <C T/2, gives [1,2] 

mo^m 

The exponential decrease of this ratio with the time distance can be traced back to 
the fact that for each gauge configuration the standard Monte Carlo allows for the 
propagation in time of all asymptotic states of the theory regardless of the quantum 
numbers of the source field O. Therefore each configuration gives a contribution to the 
signal which decreases exponentially in time, whereas it contributes 0(1) to the noise 
(variance) at any time distance. On the contrary, if in Eq. (|4.8p d is chosen large enough 
for the single thick-slice ratio to be roughly dominated by the contribution of the lightest 
state, then each factor is of order e~-^i For each configuration of the boundary fields, 
the magnitude of the product is proportional to e~^i ^, and the statistical fluctuations 
are reduced to this level. To achieve an analogous exponential gain in the computation 
of the correlation functions, the projectors T*^ have to be inserted in the proper way 
among the interpolating operators (see Ref. [9] for a more detailed discussion). 

5 Numerical simulations 

We have tested the hierarchical multi-level integration strategy described in the previous 
section for the SU(3) Yang-Mills theory by performing extensive numerical computa- 
tions. We have simulated lattices with an inverse gauge coupling of /? = 6/go = 5.7 
which corresponds to a value of the reference scale tq of about 2.93a [19,20]. The num- 
ber of lattice points in each spatial direction has been set to L = 6, 8 corresponding to 
a linear size of 1.0 and 1.4 fm respectively. For each spatial volume we have considered 
several time extents T, the full list is reported in Table [1] together with the number 
of configurations generated and the details of the multi-level simulation algorithm used 
for each run. The lattices have been chosen to test the strategy in a realistic situation 
with the computational resources at our disposal, i.e. a machine equivalent to approxi- 
matively 6 dual processor quad-core PC nodes of the last generation running for a few 
months. 
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Lattice 


L 


T 


A^'conf 




d 


Ai 


6 


4 


50 


2 


4 


A2 




5 


50 


2 


5 


A 

A3 




D 


50 


z 


D 


A4 




8 


175 


2 


4 


A5 




10 


50 


2 


5 


Ae 




12 


90 


2 


6 


A 




ID 










As 




20 


48 


3 


{5,10} 


Bi 


8 


4 


20 


2 


4 


B2 




5 


25 


2 


5 


B3 




6 


75 


2 


3 


B4 




8 


48 


2 


4 



Table 1: Simulation parameters: A'^conf is the number of configurations of the uppermost 
level, A^i ev is the number of levels and d is the thickness of the thick time-slice used for 
the various levels. 

5.1 Algorithm implementation and tests 

The basic Monte Carlo update of each link variable is a combination of heatbath and 
over-relaxation updates which implements the Cabibbo-Marinari scheme [21]. Depend- 
ing on the value of the coupling constant associated to the link at a given stage of 
the simulation, the heatbath updates the SU(2) sub-matrices by the Metropolis, the 
Creutz [22] or the Fabricius-Haan [23, 24] algorithm. In the uppermost level the gener- 
ation of the gauge field configurations consumes a negligible amount of computer time. 
At this level we perform many update cycles between subsequent configurations (typi- 
cally 500 iterations of 1 heatbath and L/2 over-relaxation updates of all link variables) 
so that they can be assumed to be statistically independent. On each of these config- 
urations we compute the "observables" P^^[T, 0], with the most expensive part being 
the estimate of the thick-slice ratio R[V^(,+(^, Vxq] at the lowest algorithmic level. The 
latter is computed by using the three-level algorithm described in the previous section, 
with the parameter values tuned sequentially level by level so to minimize the actual 
CPU cost for the required statistical precision. In all runs this has been set to be at 
most 30% of the expected absolute value of the deviation of R from 1, the latter being 
determined by some preliminary exploratory tests. As mentioned in section 14.21 the 
algorithm requires a thermalization step for each value of r which has been fixed, after 
several exploratory runs, to 500 sweeps of the full sub-lattice. 
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Figure 1: Left: the natural logarithm of RfV^p+d) shown as a function of r 

(statistical errors are smaller than symbols) for a typical configuration of the run B3. 
Right: the sum of the points in the interval [— r, r] is plotted as a function of r (one 
each eighth point for visual convenience). 



Apart from many consistency checks of the programs, we have verified several 
non-trivial properties of the basic ratios in Eqs. (j4.ip and ()4.7p . We have monitored the 
deviation from the equality 



R 



R 



-1 



(5.1) 



for several boundary configurations and all values of r, and it turns out to be com- 
patible with being a Gaussian statistical fluctuation. For the runs with d = T we 
have verified that, on each configuration and within the statistical error, the ratio 
T~ [Vt) ^o]/T[Vt, Vb] is always positive as predicted by the transfer matrix represen- 
tation. For d = T/2 the two thick-slice ratios in Eq. (j4.8p have to be equal. We have 
monitored the difference in a significant sample of our configurations, and it turns out 
to be compatible with a Gaussian statistical fluctuation as well. 

The natural logarithm of R[V^Q+(i, 14,,, r] is shown as a function of r in the left 
panel of Fig. [T] for a typical configuration of the run B3. As expected, its value is of 
0(1) for each value of r. Its almost perfect asymmetry under r — > — r, however, makes 
the sum of all the points a quantity of 0{1). This impressive cancellation, which 
is at work for T > 3 on both volumes, can be better appreciated in the right panel of 
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Figure 2: Monte Carlo history of the quantity P2 5[10,0] for the run A5. The cen- 
tral dashed line corresponds to the average value, while the other two delimit the one 
standard deviation region. 



the same Figure, where the sum of the function in the interval [— r, r] is plotted for a 
subset of values of r. It is the deviation from the exact asymmetry which flips in sign 
under a parity transformation of one of the boundary fields, and forms the signal we 
are interested in. A similar behaviour is observed for all other configurations and runs. 

The Monte Carlo history of P2'p^2[-^'^] shown in Figure [2] for the lattice A5. 
Also for all other runs we have observed reasonable Monte Carlo histories, and there- 
fore we have computed Z^/Z and its statistical error in the standard way. The run A4 
however is much noisier than the others, with rather large fluctuations due to a few 
configurations. This could be related to the fact that d = 4 is not yet large enough, and 
sizeable contaminations from the heavier states amplify the statistical fluctuations. To 
check our statistical errors, we have also carried out a more refined analysis following 
Ref. [25]. No autocorrelations among configurations have been observed, and the errors 
are fully compatible with those of the standard analysis. 

Before describing the main numerical results of the paper we mention that, for the 
runs where m = 2 is available, we have computed the quantity on the r.h.s of Eq. (13. 9p . 
As expected, it turns out to be always compatible with zero. 
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Table 2: Numerical results for various primary observables and for M (see text). 



5.2 Simulation results 

The ratios Z^ jZ have been computed for all values of m available in each run by using 
Eq. (j4.8p . The results are collected in Table [21 where they are identified by the obvious 
notation Z^^JZ. 

On each lattice the different determinations of Z^ JZ are in good agreement, and 
the sum {Z^ / Z + Z~ / Z) is always consistent with 1. For Z^ jZ a clear statistical signal 
is obtained for m = 2 only, and the larger error at m = 1 indicates that the exponential 
reduction of the noise is working as expected. To better appreciate the efficiency of the 
method, it is useful to define the quantity 



M- = --Ln 
T 



z- 



(5.2) 



whose values are reported in Table [2j With the exception of the lattice A4, it is clear 
that 0(50) measurements are enough to obtain a precision on M~ of the order of 5% on 
both spatial volumes. Sticking to the A lattices, the comparison of the relative errors 
on M~ at T = 5, 6, 10, 12 and at T = 20 indicates that the multi-level integration 
indeed achieves an exponential reduction of the noise. The most precise determination 
of Z~ jZ at each value of T is plotted in Fig. [3l Its value decreases by more than five 
orders of magnitude over the time range spanned. The symmetry constrained Monte 
Carlo clearly allows to follow the exponential decay over many orders of magnitude, a 
fact which represents one of the main results of the paper. 
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Figure 3: The quantity Z jZ as a function of T . 



The data in Table [2] confirm the expectation that at these volumes the ratio Z~ jZ 
suffers from large finite-size effects. If we enforce the theoretical prejudice that a single 
state with multiplicity 1 dominates Z~ jZ for large T, then M~ can be interpreted as 
an effective parity-odd glueball mass, which should approach its asymptotic value from 
below. Indeed this is verified at both values of L, as shown in Fig. [Hfor the A lattices. 



6 Conclusions 



The exponential growth of the statistical error with the time separation of the sources 
is the main limiting factor for computing many correlators on the lattice by a standard 
Monte Carlo procedure. The integration scheme proposed here solves this problem by 
exploiting the symmetry properties of the underlying quantum theory, and it leads to 
an exponential reduction of the statistical error. In particular the cost of computing 
the energy of the lowest state in a given symmetry sector grows linearly with the time 
extent of the lattice. 

In extensive simulations of the SU(3) Yang-Mills theory, we have observed a def- 
inite exponential reduction of the statistical error in the computation of the relative 
contribution of the parity-odd states to the partition function. The simulations needed 
at larger volumes and finer lattice spacings to provide a theoretically solid evidence for 
the presence of a glueball state, and to precisely determine its mass are now feasible 
with the present generation of computers. 
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Figure 4: The effective mass M" as a function of T. 



Since the strategy is rather general, we expect it to be appHcable to other symme- 
tries and other field theories including those with fermions as fundamental degrees of 
freedom. In QCD, for instance, the very same problem occurs already in the computa- 
tion of rather simple quantities such as the energy of the vector meson resonance, and it 
becomes even more severe for the 7]' and baryon masses. The approach presented here 
offers a new perspective for tackling these problems on the lattice. 

The integration scheme described is yet another example of how the properties 
of the underlying quantum system, namely the parity symmetry, can be exploited to 
design more efficient exact numerical algorithms for the computation of the dynamical 
properties of the theory. 
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A Numerical computation of R 

In this Appendix we describe how the ratio R[T4(,+(i, I^q, ''^]) defined in Eq. (j4.7p . has 
been computed by a three-level algorithm. The partition function T[l/j;Q_|_rf, V^^, r] is 
rewritten as 



I D4[t/]subD[C/4] e-^t^-'-l , (A.l) 



where a second temporal link Ui{yo,y) has been added to the standard degrees of 
freedom at each point of the time-slice yo = {xo + d—1). The subscript "sub" indicates 
the integration over the standard active- link variables of the thick time-slice [xq, xq + d\ 
with the spatial components Uk{x) of the boundary fields fixed to 14 (xq, x) and Vk{xQ + 
d, x) respectively. The modified action S\U^r\ reads 

5[C/,r] =5[t/] + ^(l-2r)^ReTr{c/ofc(yo,y)-C/4fc(yo,y)} , (A.2) 

where Uokiu) is defined in Eq. (j2.2p and 

Uikiy) = U^ivo, y) Ulivo + l,-y- k) ?7l(yo, y + k) Ulivo, y) . (A.3) 
If one defines the "reweighting" observable as 

0[U, r + e/2] = SU,r+e/2]-^[U,r-e/2] ^ (^^4) 

then the ratio R[T4(,+rf, T4o) '''] can be computed as its expectation value on the ensemble 
of gauge configurations generated with the action S\U^r + e/2]. In practice the aver- 
age value of the observable O is estimated by implementing the following three-level 
algorithm: 

1. Generate a thermalized configuration with the action S\U^r + e/2] by spanning 
the sub-lattice with several sweeps of the update algorithm (see section [5T]l : 

2. Compute an estimate of (O) by averaging over ng (level 0) configuration^ gener- 
ated by keeping fixed all link variables with the exception of the links Uq and C/4 
on the time-slice yoi 

3. Repeat step 2 over rii (level 1) configurations generated by keeping fixed all links 
of the sub-system with the exception of those on the time-slice yo, and average 
over the results obtained; 



^Notice that when spatial hnks are kept fixed, the set of Uo and U4, factorize and are generated 
independently. 



16 



4. Repeat step 3 over n2 (level 2) configurations generated by updating all links of 
the sub-lattice with the action £'[[/, r-|-e/2], and average over the results obtained. 

At each level the numbers uq, ni and n2 of configurations generated are chosen to 
minimize the numerical cost required to reach the desired statistical precision. Their 
values depend on d and r. In the simulations that we have carried out they range in 
the intervals no = 12 - 50, m = 50 - 120 and n2 = 50 - 300. 
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